Initialisation
## [1] "number of cores available = 1"
#Phi[1] ; eta = valeur de fin Phi[2] = valeur du noeud Phi[3] = echelle
m <- function(t, eta, phi) (phi[,1] + eta)/(1+exp((phi[,2]-t)/phi[,3]))
#=======================================#
param <- list(sigma2 = 0.05,
rho2 = 0.1,
mu = c(5,4,0.5),
omega2 = c(0.5,0.1,0.01))
#=======================================#
t <- seq(2,6, length.out = 10) #value of times
# --- longitudinal data --- #
data <- NLME_data(G = 40, ng = 12, time = t, fct = m, param = param)
getDim(data)
## G ng n N F.
## 40 12 4800 480 3
Y <- data$obs

SAEM avec simulation par MCMC
\[\log f(Y, Z ; \theta) = \langle \Phi(\theta) ; S(Y,Z) \rangle - \psi(\theta)\]
Phi <- fct_vector(function(sigma2, rho2, mu, omega2)
c(-n/(2*sigma2), -N/(2*rho2), -G/(2*omega2), G*mu/omega2), dim = c(1,1,F.,F.) )$eval
S <- fct_vector(function(eta, phi) mean((Y - get_obs(data, eta = eta, phi = phi) )^2 ),
function(eta, ...) mean(eta^2),
function(phi, ...) apply(phi^2, 2, mean),
function(phi, ...) apply(phi , 2, mean), dim = c(1,1,F.,F.) )
Metropolis Hastings
loglik.phi <- function(x, eta, Phi)
{
id <- c(1,3,4) #indice de S_1 et S_{3,.} puis S_{4,.}
sum( Phi%a%id * S$eval(eta, x, i = id) )
}
loglik.eta <- function(x, phi, Phi)
{
id <- c(1,2)
sum( Phi%a%id * S$eval(x, phi, i = id) )
}
SAEM
Initialisation
# --- Initialisation des paramètres --- #
para <- param %>% lapply(function(x) x* runif(1, 1,1.2))
para$rho2 = 0.1 ; para$mu <- c(4,5,2) ; para$omega2 <- rep(.1,3)
# --- Initialisation des chaines MC : Z_0 ---
Z <- list(eta = list( rnorm(G*ng, 0, para$rho2) %>% matrix(ncol = 1) ),
phi = list( matrix(rnorm(F.*G, para$mu, para$omega2), nrow = F.) %>% t ) )
Etape simulation et maximisation du SEAM
sim <- function(niter, h, Phih, eta, phi)
{
M <- length(phi)
eta <- 1:M %>% lapply( function(i)
MH_High_Dim_para_future(niter, eta[[i]], sd.eta(h) , loglik.eta, phi[[i]], Phih, cores = 1))
phi <- 1:M %>% lapply( function(i)
MH_High_Dim_para_future(niter, phi[[i]], sd.phi(h), loglik.phi, eta[[i]], Phih, cores = 1 ))
list(eta = eta, phi = phi)
}
maxi <- function(S)
{
list(sigma2 = S%a%1,
rho2 = S%a%2,
mu = S%a%4,
omega2 = S%a%3 - (S%a%4)^2 )
}
|
|
sigma2
|
rho2
|
mu1
|
mu2
|
mu3
|
omega21
|
omega22
|
omega23
|
|
Oracle
|
0.0494145
|
0.0989822
|
5.135294
|
3.995671
|
0.4979819
|
0.3929559
|
0.1507572
|
0.0111813
|
|
Initialisation
|
0.0511214
|
0.1000000
|
4.000000
|
5.000000
|
2.0000000
|
0.1000000
|
0.1000000
|
0.1000000
|
niter <- 25
correction.phase <- 50
MH.iter <- function(k) ifelse(k<correction.phase, 20, 10)
sd.eta <- function(k) if(k<correction.phase) .03 else .04
sd.phi <- function(k) if(k<correction.phase) c(.028, .036, .021) else .01
res <- SAEM(niter, MH.iter, para, Phi, S$eval, Z, sim, maxi,
burnin = 150, coef.burnin = 3/4,
eps = 1e-3, verbatim = 2)
saveRDS(res, 'saem.rds')
gg <- plot(res, true.value = param)
## [1] "SAEM execution time = 02min 47sec"
Result of the SAEM-MCMC
|
|
sigma2
|
rho2
|
mu.1
|
mu.2
|
mu.3
|
omega2.1
|
omega2.2
|
omega2.3
|
|
Real value
|
0.0500
|
0.1000
|
5.0000
|
4.0000
|
0.5000
|
0.5000
|
0.1000
|
0.0100
|
|
Estimated value
|
3.3205
|
0.1738
|
3.8711
|
5.0436
|
2.0748
|
0.2983
|
0.4205
|
0.1780
|
|
Rrmse
|
65.4090
|
0.7377
|
0.2258
|
0.2609
|
3.1495
|
0.4034
|
3.2051
|
16.7954
|


